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Non-linear k-e-v 2 modeling 
with application to high-lift 

By F. S. Lien 1 AND P. A. Durbin 2 


The k-e-v 2 model has been investigated to quantify its predictive performance on 
two high-lift configurations: 2D flow over a single-element aerofoil, involving closed- 
type separation; 3D flow over a prolate spheroid, involving open-type separation. A 
‘code- friendly’ modification has been proposed which enhances the numerical stabil- 
ity, in particular, for explicit and uncoupled flow solvers. As a result of introducing 
Reynolds-number dependence into a coefficient of the e-equation, the skin-friction 
distribution for the by-pass transitional flow over a flat plate is better predicted. 
In order to improve deficiencies arising from the Boussinesq approximation, a non- 
linear stress-strain constitutive relation was adopted, in which the only one free 
constant is calibrated on the basis of DNS data, and the Reynolds- stress anisotropy 
near the wall is fairly well represented. 


1. Introduction 

Eddy-viscosity models based on the linear Boussinesq relations are known to be 
afflicted by numerous weaknesses, including an inability to capture normal stress 
anisotropy, insufficient sensitivity to secondary strains, seriously excessive genera- 
tion of turbulence at impingement zones, and a violation of realizability at large 
rates of strain. Notwithstanding these defects, eddy-viscosity models remain popu- 
lar, and their use in complex flows is widespread due, principally, to their formalistic 
simplicity, numerical robustness, and computational economy. Second-moment clo- 
sure, on the other hand, accounts for several of the key features of turbulence that 
are misrepresented by linear eddy- viscosity models, but is considerably more com- 
plex and can suffer from poor numerical stability due to the lack of dominance of 
second-order fragments in the set of terms representing diffusion. As a result, the 
CPU requirements for second-moment closure models can be high, especially in 3D 
flows. 

A potential alternative to second-moment closure, but one which retains advan- 
tageous elements of the linear eddy-viscosity framework, is to use a constitutive 
relation that equates the Reynolds-stresses to a non-linear expansion in powers of 
the mean rate of strain and rate of rotation tensors. This may be cast in the form 
of a sum of terms, each pre-multiplied by an apparent viscosity — hence the term 
‘non-linear eddy-viscosity models’. Examples include the models of Speziale (1987), 
Shih et al. (1993), Durbin (1995a), Craft et al. (1995) and Lien et al. (1996). The 
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main differences between the above modeling strategies can be summarized in the 
following table: 


Authors(s) 

Model 

form 

Order in the stress- 
strain relationship 

Number of turbulence 
transport equations 

Speziale 

(1987) 

High-Re 

quadratic 


Shih et ai 

(1993) 


quadratic 

2, k-e 

Durbin 

(1995a) 


KBS9H1 

3, k — e — v 2 


m 


cubic 

3 , k - e — A 2 


(1996) 

mssam 


2, k - £ 


The ,42 value — the second Reynolds-stress invariant — in Craft et a/.’s k — e — A 2 
model is obtained by solving a related transport equation as follows: 

0, A 2 + V • V.4 2 = -2^-(d k + P k - e) 
k 

+2-±{d ij + P l) + <f>ij — £ij), (1) 

with fragments consistent with second- moment closure. In order to be free from 
topological constraints, the unit vector in the wall- reflection term is replaced by the 
length-scale gradient. The expansion of (1) in 3D curvilinear coordinate systems is 
tedious and prone to error. Also, a major drawback of this model is the high level 
of sensitivity to the near-wall grid parameters, including resolution, distribution, 
and aspect ratio. 

The u 2 -equation in Durbin's k — e — v 2 model, to be addressed in Section 2, was 
simplified from second-moment closure on the basis of the IP pressure-strain model 
in conjunction with elliptic relaxation. This approach is algorithmically simple, 
applicable to the low-Re region, and naturally mimics the kinematic blocking effect 
on the turbulence of a solid wall. 

Another important feature which distinguishes Durbin’s model from most others 
is the expression of eddy- viscosity v u which plays an important role in determining 
the correct level of shear stress. In Craft et aV s model, 


v t = 0.734- 


,[1 - exp {—0.145 exp(1.3yy 5 ^ 6 )}] - 0.8exp(-.R t /30)} 


1 + 1.8 T) 


1 + 0.6.4 2 + 0.2AI 5 


(kT) (2) 


where 


r t) — 1 T 


1 — exp(- 


4 3 

_ 2 _ > 

0.125' 


1 + 


4yJexp(-R t /20) 


, r\ — max(5, (3) 


and S and Q are strain and vorticity invariants. While in Durbin’s model, 


= 0.19— (fcT). 


(4) 
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One distinct difference between Eqs. (2) and (4) is that the latter does not require 
any damping function: a result of using v 2 as the velocity scale in the direction of 
the wall. The former, on the other hand, sensitizes v t to S, ft, R t (i.e., Reynolds 
number) and A 2 , with the functional dependency being carefully calibrated on a 
range of flows, including straining flow, channel flow, impinging jet, and transi- 
tional flow. However when this model was tested for turbomachinery flows at (and 
near) off- design conditions, the size of the leading-edge separation bubble was over- 
estimated, and in some cases no converged solution could be obtained. This is due 
to 7 ) (strain and vorticity) and A 2 being too large along the curved shear layer . As 
a result of both parameters appearing in the denominator of v t expression, the level 
of shear stress was significantly under-predicted (Chen, 1996). 

In the present work, the k-e — v 2 model of Durbin (1995b) is applied to high-lift 
configurations, both 2D and 3D. In the course of this study, numerical instability 
arising from the boundary condition atjvall was encountered, due to our use of a 
solution algorithm that uncouples the v 2 and /-equations. A ‘code- friendly’ mod- 
ification is introduced, which not only circumvents this numerical difficulty, but 
also gives better predictions for transitional flows. This variant is then combined 
with the non-linear stress-strain constitutive equation with the aim of improving 
the near-wall behavior of normal-stress anisotropy. 

2. k — £ — v 2 model 

The turbulence model uses the standard k — e equations: 

dtk + U'Vk = Pk-e + [(v+—)Vk], (5) 

<7k 

d t s + U-Ve= C( ' Pk ~ — + [(i/ + — )Ve]. (6) 

1 <T e 

On no-slip boundaries, y — > 0, 

k = 0, e — > 2v^t. (7) 

r 

The v 2 transport equation is 

d t v 2 + U • Vt? 2 = kf — nv 2 j + V • [{v + ^)Vr 2 ], (8) 

ft 


where kf represents redistribution of turbulence energy from the streamwise com- 
ponent. Non- locality is represented by solving an elliptic relaxation equation for 
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The Boussinesq approximation is used for the stress-strain relation: 

_ TiJuJ 2 _ v t 

a,,-— 3 k 

where the eddy viscosity is given by 

v t = cyr. 


(ii) 


(12) 


The constants of the model are: 


C ^ = 0.19, <Jk = 1, cr € = 1.3, 

C ( i = 1.55, C< 2 = 1-9 

Ci = 1.4, C 2 = 0.3, C L = 0.3, C v = 70. 

As y — ► 0 — y being the minimum distance to walls — and k 
Eq. (8) becomes: 

vdlv 2 —2 nv%r = kf. 

y yl 


(13) 
0 \j2v)ey l , 


(14) 


The viscous and kinematic conditions at the wall show that v 2 should be 0(y 4 ) as 
y -y 0. In the original k — e — v 2 model, n = 1, yielding the boundary condition for 


/ 


x/m (24 — 4n)v 2 v 2 , _ 20 v 2 v 2 

} " e(o) y 4 |n=1 = ~7(0W ' 


(15) 


on no-slip walls. 


2.1 Code-friendly modification 

Equation (15) works fairly well for coupled, implicit solvers [e.g, INS2D code 
of Rogers &: Kwak (1990)]. However, for explicit and uncoupled schemes, numer- 
ical instability arising from y 4 in the denominator of Eq. (15) sometimes occurs. 
Therefore, a code-friendly modification is made here by setting n = 6, which allows 
/( 0) = 0 to be imposed as the boundary condition. In addition, C t i and C € 2 are 
replaced by 

C el = 1.55 + exp(-A f ^)U <=0 .00285, C <2 = 1.92, (16) 


where R y = yy/k/u , and the other model constants 


are: 


C^ = 0.19, a* = 1, o t — 1.5, 

Ci = 1.4, C 2 = 0.3, C L = 0.17 C, = 70. (17) 

Fully- developed channel flow 

The model constants, in particular A t = 0.00285 and Cl = 0.17, were first 
calibrated with the channel-flow DNS data of Kim e< al. (1987) and then optimized 
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FIGURE 1. Channel flow: (Left) mean velocity; (Right) k and v 2 . DNS: • velocity; 
A k; □ v 2 . 

on the basis of 2D/3D separated flows to be presented later. As seen in_Fig. 1, 
both the mean-velocity and turbulence profiles, the latter including k and v 2 , agree 
reasonably well with the data. 

2.1.2 By-pass transitional flow over a flat plate 

The second case examined here is the flow over a flat plate with free-stream turbu- 
lence intensity T u = 3% and dissipation length scale = 10 mm. The experimental 
study was conducted at Rolls Royce Aeroengines in Derby, UK. The skin-friction 
distributions, obtained with the original and code-friendly k - e — v 2 variants and 
Launder- Sharma model (1974), are shown in Fig. 2(L). As seen, introducing the 
R y -dependency in C € \ for the code-friendly variant improves transition predictions. 
Although the resulting onset of transition is slightly earlier than that returned by 
the Launder- Sharma model, the length of transition is better represented. As the 
flow becomes fully turbulent, the velocity profiles obtained with both k — e — v 2 
variants are almost identical as demonstrated in Fig. 2(R). 

2.2 Non-linear constitutive relation 

A general constitutive relation of the type proposed by Pope (1975) can be written 
as: 

0 10 

ay = U ~t L ~ i Sij = YG\S ij ,Sl i] ,vyk,T)T t ). (18) 

k 6 \=i 

where = S i} , T? = S ik £l k j - tt, k S k} , T? } = S lk S k j - ^ tJ Si k S k i ■ ■ ■. Truncating 
at the third term for simplicity gives rise to 

ai j = -jS tJ + G 2 (S tk n k} - Q ik S k] ) + G\S tk S k] - l*yS,*S«). (19) 


, _dU 1 dUj_ dU, dUj 

,J dx j dxi ' dxj dx{ 


where 


( 20 ) 
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log, (C r ) 



4.0 4.5 5.0 5.5 6.0 -1 0 1 2 


FIGURE 2. Flat plate: (Left) skin friction; (Right) mean-velocity profile. • expt.; 
k — e\ o original k — e — v 2 ; modified k — e — v 2 

Two constrains for parallel flow will be imposed: 


«22 — an = OLCI 22 = oca- 7, 

where a—f — ^ These yield 


G 2 



- (1 + a ) a y2 rp 2 

2 S 2 


( 21 ) 

( 22 ) 


where 5 = j|-~| or (= jy/SijSij :/2, in general) and T is defined in Eq. (6). The 
remaining unknown, a, can be evaluated from DNS data of channel flow (Kim et 
al . , 1987), boundary-layer flow (Spalart, 1988) and flow over a backward-facing step 
(Le et a/., 1993). As seen in Fig. 3, 



fits DNS data reasonably well. The algebraic model was initially used by Durbin 
(1995a) as an a postiori formula for evaluating ujuj. In order to apply Eq. (19) to 
mean flow prediction while preventing computational intractability, the coefficients 
G 2 and G 3 are modified as: 


G 2 


1 (1 a ) a v 2 rp2 _ 3 (1 + a ) a v 2 2 

4 S 2 + 1 ’ “ 2 S 2 + 1 


(24) 


3. Numerical method 

All flows have been computed with the STREAM general geometry, block-struct- 
ured, finite-volume code (Lien k Leschziner, 1994a). Advection is approximated 
by a TVD scheme with the UMIST limiter (Lien k Leschziner, 1994b). To avoid 
checkerboard oscillations within the co-located storage arrangement, the ”Rhih and 
Chow 1 interpolation method (1983) is used. The solution is effected by an iterative 
pressure-correction SIMPLE algorithm, applicable to both subsonic and transonic 
conditions (Lien k Leschziner, 1993). 
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FIGURE 4. A-aerofoil: geometry and partial grid 

4. Results and discussion 

4-1 Aerospatiale A-aerofoil 

Computations for the A-aerofoil have been performed at 13.3° incidence, with 
transition on the suction side prescribed at 12% of chord. The geometry and a 
partial view of the grid are given in Fig. 4. The Reynolds number, based on free- 
stream velocity and chord length, is 2.1 x 10 6 . Solutions have been obtained on a 
grid containing 177 x 65 lines, extending to 10 chords into the free stream. 

In total, four turbulence-model variants have been applied to this case [compar- 
isons to second-moment closure can be found in Lien & Leschziner (1995)]: 

(1) the low-Re k — e model of Lien Sz Leschziner (1993); 

(2) the original k — £ — v 2 model of Durbin (1995b); 

(3) the code-friendly variant; 

(4) the above variant combined with the non-linear stress-strain relation. 

The skin-friction and wall-pressure distributions obtained with three linear eddy- 
viscosity models, one k — e and two k — e — v 2 , are compared in Fig. 5. These, as 
well as the associated profiles of streamwise velocity and shear stress on the suction 
side in Figs. 6-7, clearly demonstrate the superiority of k — e — v 2 variants relative 
to the conventional k — e model. 

Attention is turned next to comparisons between linear and non-linear k — e — v 2 
models in Figs. 8-10 for profiles of streamwise velocity and Reynolds normal- 
stresses. It is found from these figures that the Reynolds-stress anisotropy is fairly 
well predicted by the non-linear model at x/c=0.5, which is consistent with the con- 
straints in Eq. (21) imposed on the constitutive equation. As the flow approaches 
the trailing edge, streamline curvature arising from secondary strain becomes impor- 
tant and the omission of its production term (~ in the v 2 -equation is no longer 
valid, resulting in large discrepancies between predictions and data at x/c=0.9. 

4-2 DLR prolate spheroid 

The shape of this body and a partial view of the numerical grid surrounding it 
are shown in Fig. 11. The Reynolds number, based on the chord, is 6.5 x 10 6 . 
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FIGURE 5. A-aerofoil: (Left) skin friction; (Right) pressure coefficient. • expt.; 
k — e; original k — e — v 2 ; modified k — e — v 2 



Computations have been performed at 30° incidence in which transition is free. 
The solution domain, containing 65 x 65 x 65 lines, extends 10 chords into the outer 
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y/c x/c= 0.83 y/c x/c= 0.90 y/c x/c= 0.96 



0.0 0.005 0.0 0.005 0.0 0.005 

y/c x/c= 0.30 y/c x/c= 0.50 y/c x/c= 0.70 



0.0 0.005 0.0 0.005 0.0 0.005 


FIGURE 7. A -aerofoil: profiles of shear stress, o expt.; k — e; original 

k — e — v 2 ; modified k — e — v 2 

stream. 

Numerical solutions have been obtained with two models [comparisons with second- 
moment closure can be found in Lien &; Leschziner (1995b)]: 

(1) the low-Re k — e model of Lien &: Leschziner (1993); 

(2) the code-friendly k — e — v 2 variant in conjunction with Launder and Kato’s 
modification in the turbulence production P* (1993). 

A well-known defect of any conventional, linear eddy- viscosity model is that it 
predicts excessive levels of turbulence energy in impingement regions, due to the fact 
that the irrotational strains appearing in the turbulence-energy equation (~ S^Sjj) 
act to generate turbulence irrespective of their sign. The rationale behind Launder 
& Kato’s proposal is to partially replace the strain by the vorticity, i.e. 

P k =0.5i/,Syn y . (25) 

A similar idea, based on ^realizability’ constraints on the turbulence time scale, 
has been sugge sted rece ntly by Durbin (1996), in which a upper bound to k/e 
proportion to yj2j S lJ S lJ was introduced. As a result, the rate of turbulence-energy 
generation in the vicinity of stagnation regions becomes linear , which is similar to 
that returned by most of the non-linear eddy- viscosity models mentioned in Section 
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FIGURE 8. A-aerofoil: profiles of streamwise velocity, o expt.; linear 

k — e — v 2 ; non-linear k — e — v 2 



0.0 0.025 0.0 0.025 0.0 0.025 


FIGURE 9. A-aerofoil: profiles of streamwise normal stress, o expt.; linear 

k — e — v 2 ; non-linear k — e — v 2 



FIGURE 10. A-aerofoil: profiles of transverse normal stress, o expt.; linear 

k — e — v 2 ; non-linear k — e — v 2 
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FIGURE 11. Prolate spheroid: geometry and partial grid 


1; those models use 


C\ 




o 

l + .dv^Sp"' 


(26) 


Azimuthal variations of skin-friction magnitude and direction at four streamwise 
locations are shown in Figs. 12 and 13, and the circumferential distributions of wall 
pressure are given in Fig. 14. As seen at x/2 a ~ 0.223. the k — s — r 2 model in 
conjunction with Launder & Kato’s modification returns a transition-like behavior 
in the boundary layer close to the windward side. Although the model is unable 
on fundamental grounds to predict any aspect of natural transition, the predicted 
transitional phenomenon is mainly due to a strong suppression of turbulence energy 
at the impingement regions, in which the flow becomes ‘laminar \ combined with the 
fact that the free-stream turbulence diffuses into the boundary layer and ultimately 
triggers transition. It is clear from Fig. 14 that the extent of pressure 1 plateau 
regions, signifying the azimuthal extent of separation zone, at x/2 a > 0.564 are 
under-estimated by both models. This observation is consistent with the azimuthal 
distributions of skin-friction direction 7 shown in Fig. 13; 7 — 0 denotes either 
the separation or the reattachment point. The performance of k — s — v 2 model 
is slightly better than that of k — e in terms of the extent of the separation zone. 
Some of the discrepancies between predictions and experiment might be due to the 
grid density adopted here; in particular, close to the rear end of the spheroid it is 
too coarse and a grid-refinement test is required. 


5. Conclusions 

A computational study has been undertaken to investigate the predictive capa- 
bilities of k — e — v 2 variants when applied to high-lift configurations, including 
2D aerofoil and 3D prolate spheroid. Both the linear and non-linear stress-strain 
constitutive relations are examined. The outcome of the present study may be 
summarized as follows: 
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modified k — e — v 2 ; 


(1) The k — e-v 2 model and its variants, whether linear or non-linear, return superior 
predictions relative to the conventional k — e model. 

(2) This superiority can be attributed to the use of v 2 as the velocity scale in the 
eddy-viscosity expression without resorting to an ad hoc damping function. 

(3) The v 2 is obtained from a simplified form of Reynolds-stress transport equation, 
governing the turbulence intensity normal to streamlines, the pressure-strain term 
of which is represented mathematically by an elliptic relaxation model. 

(4) A code-friendly modification is proposed here, including the assurance of the 
near- wall behavior v 2 — > 0(y 4 5 6 ) as y -> 0, the introduction of R y -dependency in 
C € i, and the use of / = 0 as the boundary condition on no-slip boundaries. As a 
result, the numerical stability, in particular, for the uncoupled solution procedure 
used herein is greatly enhanced. 

(5) The introduction of R y in C t \ yields improved results for the transitional flow. 
However, it requires the minimum distance to walls, which can be difficult to 
apply to complex geometries. 

(6) Following a similar idea suggested by Durbin Sz Laurence (1996), a first attempt 
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FIGURE 13. Prolate spheroid: skin- friction direction, o expt.; 

modified k — e — i> 2 ; 


has been made by adopting 

C ( i = 1.44(1 + 0.0333 yfkft?). 


C t 2 = 1.85, C L = 0.188, 


and preliminary results for flows over a flat plate and the A-aerofoil, described in 

Sections 2. 1.1-2. 1.2, are given in Figs. 15-17. As seen, the use of yj k/v 2 returns 
very similar mean-velocity profiles for the A-aerofoil case. However, the onset of 
transition for the flat-plate case is too early and the length of transition is too 
long. _ 

(7) In order to improve the performance of k — e — v 2 model for both transitional and 
fully turbulent flows, in particular, in complex geometries, instead of adopting 

R y and yj k/v 2 , there is a need to devise a new parameter, depending on the local 
Reynolds number and avoiding the use of the minimum distance to walls. 

(8) The level of normal stress anisotropy returned by the non-linear model is fairly 

well represented at the mid-chord of A-aerofoil, where the curvature effect is 
unimportant. Close to the trailing edge, however, both t£ 2 , v 2 and, consequently, 
k and its production P* are under-predicted. Since Pk = + * * * 

and the mean-velocity profile and, hence, its gradient at x/c = 0.9 are in good 
agreement with the data, this indicates that v t is too low, which is consistent 
with the under estimation of v 2 at the same location. 
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-i J 

FIGURE 14. Prolate spheroid: pressure coefficient, o expt.; □ inviscid solution; 
k — e; modified k — e — v 2 ; 



(9) For open 3D separation, the size of separation zone, reflected by the_ azimuthal ex- 
tent of pressure plateau, is slightly under-predicted by the k — e — v 2 model, which 
might be partially attributed to the grid density adopted here being insufficient. 
(10) To ensure a wide range of applicability of the non-linear model, the free coefficients 



FIGURE 16. A-aerofoil: profiles of streamwise velocity, o expt.; 
R y ; based on \J k/v 2 


based on 


lOOOxCf 



lOOOxCf 



FIGURE 17. Prolate spheroid: skin-friction magnitude, o expt.; based on 

R y ; based on y k/v 2 

and their associated functionals need to be more carefully optimized by reference 
to different types of flow, featuring separation, impingement, swirl, rotation, and 
transition. 
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